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BEAM-PLASMA  HEATING  MODEL  CODE  AND  COMMENTARY 


INTRODUCTION 

This  code  listing  and  commentary  is  intended  to  supplement  the  paper, 
"Beam-Plasma  Heating  Model,"  by  D.  A.  Hammer,  A,  W.  All  and  myself.1 
Briefly,  the  paper  describes  a one-dimensional  ionization  and  heating 
model  as  it  applies  to  results  of  several  electron  beam-plasma  inter- 
action experiments.  Beam  energy  is  deposited  resistive ly  in  the  plasma 
at  a rate  T\)a,  where  J is  the  return  current  density  and  T]  the  plasma 
resistivity,  both  classical  and  anomalous,  due  to  ion-acoustic  or 
electron-electron  mode  turbulence.  Principle  energy  losses  include 
ionization,  line  radiation,  inelastic  electron  impact  excitation, 
bremmst rah lung  and  radiative  recombination.  The  level  of  ionization  and 
plasma  heating  are  computed  as  a function  of  neutral  gas  pressure, 
beam  rise-time,  pulse  width  and  current  density,  and  resistivity  model. 
Plasma  dynamic  and  kinetic  effects  such  as  expansion  and  end  loss  are 
not  explicitly  included  in  the  model. 

The  computational  model  assumes  that  the  avalanche  breakdown  pro- 
cess driven  by  the  large  electrostatic  and  inductive  electric  fields  at 
the  beam  front  has  been  completed.  Previous  theoretical  studies"  have 
determined  that  the  resultant  plasma  has  a density  of  the  order  of 

1014  cm-3  and  a temperature  of  a few  eV.  Plasma  conductivity  is  then 
Not*:  Manuscript  lubmittcd  March  5,  1979. 
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Sufficiently  high  that  neutralization  of  the  beam  front  Inductive 
electric  field  drives  a backstreaming  plasma  "return  current"  vlthln 
the  beam  channel.3  When  the  beam  Is  present,  the  plasma  current  density 


J is  given  by  J “ — ).,  where  J.  Is  the  beam  current  density.  When  the 

p p D D 

beam  pulse  Is  over  J a 0. 

P 

The  remainder  of  this  paper  consists  of  a listing,  a definition  of 
terms,  and  a CSN  step  cooanentary  of  the  beam-plasma  heating  code  set 
forth  in  Section  II  of  Ref.  1.  This  report,  along  with  Ref.  1,  should  be 
sufficient  to  allow  further  application  of  this  heating  code  to  a wide 


and  varied  range  of  beam  and  plasma  parameters  and  resistivity  models. 


CODE  LISTING 

This  code  calculates  the  density  and  temperature  in  a "return  cur- 
rent" heated  plasma  in  the  presence  of  a large  initial  neutral  fraction. 
Avalanching  is  assumed  completed  when  the  calculation  starts  (i.e.,  no 
electric  field  is  applied).  The  energy  equations  (Eqs.  (9)  and  (11), 

Ref.  1)  calculaced  in  subroutine  FINI,  ENTRY  F,  Include  Spitzer  and 
electron  neutral  resistivities  and  three-body  recombination  as  heat 
souses.  Loss  mechanisms  are  line  radiation  and  ionization.  The 
density  equation  (Eq.  (10),  Ref.  1)  also  calculated  in  ENTRY  F,  Includes 
collisional  ionization  and  two  and  three-body  recombination.  Integration 
of  the  energy  and  density  equations  is  carried  out  by  subroutine  LNT 
(Boris  and  Windsor).4 
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DEFINITION  OF  TERMS 


■ 


ANOME 

ANOMEE 

ANOMI 


Electron  portion  of  resistivity  as  expressed 
in  Eq.  (12a),  Section  II,  of  Ref.  1. 

Electron-electron  resistivity  as  expressed  in 
Eq.  (4). 

Ion  portion  of  resistivity  as  expressed  in 
Eq.  (12b) . 


ANOMIA 

ANOMLH 

ANOMWV 
BKG 
8 (KG) 

BREMS 

CLAS  RES 

CLASSR 

CL  ION  HT 

COULOG 
CRDENT 
CRD  NT 


ion-acoustic  resistivity  as  expressed  in  Eq.  (3) . 

Lower  hybrid  resistivity  not  specifically  given 
in  Ref.  1,  but  operative  in  some  of  the  runs. 

Wave  resistivity  as  expressed  in  Eq.  (8)  . 

Applied,  axial,  magnetic  guide  field  in  kG. 

Format  designating  printout  of  initial  value  of 
magnetic  field  BKG. 

P„  in  Ref.  1,  the  plasma  power  density  in 
ev/cm3  sec  lost  to  bremsstrahlung. 

Format  designating  printout  of  instantaneous 
classical  resistivity  value. 

Classical  resistivity,  as  expressed  in  Eqs. 

(la)  , (lb)  and  (2) . 

Format  designating  printout  of  power  density 
used  for  classical  ion  heating. 

Coulomb  logarithm,  ln/\  • 

Instantaneous  beam  current  density  in  kA/cm2. 

Format  designating  printout  of  instantaneous 
beam  current  density  in  kA/cm2. 


CS 


Ion  sound  velocity. 


CUM  ENERGY 


CURDEN 


Format  designating  printout  of  total  accumulated 
energy  density  transferred  from  the  beam  to  the 
plasma  system,  in  eV/cm. 3 

Assymptotic  beam  current  density  in  kA/cm2. 


» 


i 


3 


"■j*. 


CUTOFF 


DENE 


Ratio  of  instantaneous  beam  current  density 

(CRDENT)  and  a minimum  current  density  (STABU) 

or  J . as  expressed  in  Eq.  (6) . 
min 

Instantaneous  value  of  the  electron  density 


DENEUT 


Instant aneous  value  of  the  neutral  density 
cm*3. 


DENI 

DENMAX 

DENS ITY 


Instantaneous  value  of  the  ion  density,  equal 
to  DENE  for  a hydrogen  plasma. 

Neutral  hydrogen  fill  density  cm*3.  Equal  to 
DENE  plus  DENEUT. 

Format  designating  printout  of  initial  value 
of  electron  density  (Y20)  in  cm”3. 


DY(1) 


Electron  energy  equation  as  expressed  in  Eq.  (9) • 


DY(2) 

DY(J) 

EDENS 

EE  RES 

ENERGY 

HO, HOI 


Electron  (and  ion)  density  equation  as  expressed 
in  Eq.  (10). 

Ion  energy  equation  as  expressed  in  Eq.  (11) . 

Instantaneous  value  of  electron  density  times 

1 x 10*15. 

Format  designating  printout  of  electron-electron 
resistivity. 

Total  beam  energy  density  in  eV/cm3  deposited 
resistively  in  the  plasma  at  a rate  T\j2. 

Calculation  time  step  length. 


I AC  RES 

IH01 

IMAX 

IONIZ  RAD 


Format  designating  printout  of  ion  acoustic 
resistivity. 

Time  step  length  times  1 x 10®. 

The  number  of  time  steps  printed  out. 

Format  designating  printout  of  beam  power 
density  in  eV/'crn3  sec  used  for  ionization. 


ITIME  Time  step  length. 

JMAX  The  number  of  different  current  densities  to 

be  calculated. 
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LINE  RAD 


Format  designating  printout  of  the  plasma 
power  density  in  eV/cm3  sec  lost  through  line 
radiation. 


LWHYB  RES 

Format  designating  printout  of  the  lower  hybrid 
resistivity. 

MAX  J-BEAM 

Format  designating  printout  of  initial  value 
of  the  peak  beam  current  density  in  kA/cm2. 

NEUTRAL  DENSITY 

Format  designating  printout  of  initial  value 
of  the  neutral  hydrogen  fill  density  cm'3. 

QI 

The  classical  collisional  electron-ion  energy 
transfer  rate. 

RAD ION 

Beam  power  density  in  eV/cm3  sec  used  for  ion- 
ization. 

RADLIN 

Plasma  power  density  in  eV/cm3  sec  lost  through 
line  radiation. 

RATIO 

Ratio  of  plasma  electron  drift  velocity  VD  to 
the  critical  velocity  VCRIT.  RATIO  is  a 
measure  of  whether  the  ion  acoustic  instability 
is  expected.  Yes,  i ' RATIO  > 1. 

RESIST 

Sum  of  all  the  resistivities  employed  in  the 
model . 

S 

Ionization  rate  coefficient  as  expressed  in 

Eq.  (13). 

SI 

Excitation  rate  coefficient  as  expressed  in 

Eq.  (HO. 

STABLJ 

Minimum  current  density,  j . as  expressed  in 
Eq.  (6).  mln 

T,T(NS) 

Format  designating  printout  of  time  in  nsec  of 
the  calculation. 

TEMPERATURE 

Format  designating  printout  of  initial  value  of 
the  electron  temperature  in  eV. 

TE 

Format  designating  printout  of  electron  temper- 
ature in  eV. 

TEO 

Initial  electron  temperature  in  eV. 
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THREE B 

Three-body  recombination  rate  coefficient, 
orf  as  expressed  in  Eq.  (.16), 

TI 

Format  designating  printout  of  ion  temperature 
in  eV. 

TIME  STEP 

Format  designating  printout  of  initial  value 
of  time  step  length  (HO)  in  nsec. 

T PULSE 

Beam  pulse  width  in  nsec. 

TWOBOD 

Two-body  recombination  rate  coefficient,  a-,  as 
expressed  in  Eq.  (15). 

VCRIT 

Critical  velocity  as  expressed  in  Eq.  (5).  (See 
RATIO) . 

VD 

Plasma  electron  drift  velocity. 

VE 

Plasma  electron  thermal  velocity. 

VTHI 

Ion  thermal  velocity. 

Y(l) , 3/2NKT 

Defined  as  ?/2  the  initial  electron  density 
(Y20)  times  the  initial  electron  temperature 
(TEO) , in  eV/cm3. 

Y 20,Y(2) 

Initial  electron  density,  cm'3. 
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CSN  STEP  COMMENTARY 

The  following  commentary  la  intended  to  point  out  specific  areas 
of  importance  in  understanding  and  utilizing  this  computational 
program.  References  are  made  to  CSN  numbers  listed  down  the  left 
hand  side  of  the  code  listing. 

A.  Main  Program: 

CSN  0007  Initial  data  is  read  in  at  this  statement. 

The  present  number  of  time  steps  (IMAX'  is  100. 

The  total  number  of  beam  current  densities 
(JMAX)  is  presently  6,  (2, 3*1*  ,b,  10  and  1*5  kA/cm^. 
The  initial  electron  temperature  (TK01  is  taken 
as  2.‘;  eV,  and  is  determined  beforehand.  The 
Initial  electron  density  (Y 201  is  also  estimated 
beforehand,  and  is  usually  betveen  0.1  and  l'1<  of 
DENMAX,  the  neutral  hydrogen  fill  density. 

Present  values  of  Y20  range  between  X \ lO^cm"'’ 
and  3 x 1014cm"'3.  Values  of  DENMAX  range  between 

0. 7  x 1015cm*1  (lOOmT)  and  3.3  x 10u'cm*  ' (s00mT> . 
The  calculation  time  step  length  (HO)  is  presently 
2 x 10"9  sec.  The  first  50  steps  are  taken  at 
SH0  or  10  nsec  Increments.  See  CSN  00$3  through 
00°'},  especially  OO'Jk  and  00^8,  for  setting  the 
time  increments.  The  time  of  calculation  TINS! 

is  printed  out  down  the  left  hand  side  of  the 
computation  printout,  see  CSN  0027  for  this  for- 
mat statement.  The  beam  pulse  width  (TPULSE1  is 
presently  set  at  70  nsec  for  our  accelerator. 

The  rise  and  fall  times  are  set  in  CSN  00X2-0036. 
They  are  presently  10  nsec.  It  is  also  possible 
to  employ  a step  function  beam  pulse  by  specify- 
ing Just  the  rise  time.  In  that  case.  TTVLSE  is 
disregarded.  The  applied,  axial,  magnetic  guide 
field  (BKG)  is  set  at  2 kG. 

CSN  0020  The  assymptotic  beam  current  densities  are  read 

in  at  this  statement.  One  data  card  is  required 
for  each  value  of  Cl’RDEN. 

CSN  0023  This  is  the  format  statement  for  the  display  of 

the  initial  values  of  the  run  as  they  are  listed 
across  the  top  of  each  section  of  the  computation 
printout.  These  Include  TIME  STEP,  or  HO  in  nsec, 

1 .  NET  or  Y(  1)  in  eV/cm\  DENSITY,  the  initial 
electron  density  Y20  or  Y(2)  in  cm*  \ TEMTERATVRE . 
the  initial  electron  temperature  TE0  in  eV,  MAX 
J-BEAM,  the  initial  value  of  the  peak  beam  current 


density  CURDEN  in  kA/cm2,  NEUTRAL  DENSITY,  the 
initial  value  of  the  neutral  hydrogen  fill  density 
DENMAX,  in  cm'3,  and  B(KG)  the  initial  value  of 
the  magnetic  guide  field  BKG  in  kG. 

CSN  0027  This  is  the  format  statement  for  the  headings 

of  the  columns  of  computational  results  listed 
in  the  printout.  From  left  to  right  we  have: 

T(NS)  in  nsec,  the  time  from  the  start  of  the 
beam  pulse  to  the  time  of  the  calculation,  the 
electron  density,  DENSITY,  in  units  of  10l5cra'3, 
the  electron  and  ion  temperatures,  TE  and  TI, 
respectively,  the  Instantaneous  beam  current 
density,  CRDNT,  in  kA/cm?,  RATIO,  the  ratio 
of  plasma  electron  drift  velocity,  VD,  to  the 
critical  velocity,  VCRIT,  as  defined  in  Eq.  (t;) 
of  Ref.  1.  RATIO  is  a measure  of  whether  the 
ion  acoustic  instability  is  expected.  If  RATIO 
> 1,  the  ion  acoustic  mode  is  operative.  IONIZ 
RAD  is  the  beam  power  density  in  eV/cm3sec  used 
for  ionization.  LINE  RAD  is  the  plasma  power 
density  lost  through  line  radiation.  CL  ION  HT 
is  the  power  density  used  for  classical  ion 
heating.  IAC  RES,  LWHY3  RES,  CLAS  RES,  and  EE  RES 
in  the  last  column,  are  the  ion  acoustic, 
lower  hybrid  classical  and  electron-electron 
resistivities  respectively  as  described  in 
Section  II  of  Ref.  1.  They  have  the  units  of 
(eV/cm3sec) / (kA/cm2)  2,  so  that  power  density  is 
in  eV/cm3  sac  when  current  density  is  in  kA/cm2. 
BREMS  is  the  plasma  power  density,  P in  Ref.  1, 
lost  to  bremsstrahlung.  CUM  ENERGY  is  the  total 
accumulated  energy  density  transferred  from  the 
beam  to  the  plasma  system,  in  eV/em3. 

CSN  0028  This  statement  calls  for  subroutine  FINI  which 

calculates  the  energy  equations  given  in  Eqs.  (^) 
and  (11)  and  the  density  equation  given  in  Eq. 

(10)  in  Ref.  1.  The  common  variables  are  listed 
within  the  parenthesis. 

CSN  0C?2-003b  These  statements  specify  the  beam  pulse  rise 

(0035)  and  fall  (0033)  time,  in  this  case  being 
10  nsec.  To  specify  a step  function  pulse,  only 
statement  0035  would  be  necessary. 

CSN  0038-0080  Here  begins  the  calculation  of  quantities  for 
the  printout.  All  of  these  quantities  are  de- 
fined in  Section  III  of  this  paper.  Note  the 


13 


ord«r  of  CSN  0054  and  00*:5.  Sine*  the  computer 
"remembers"  the  lest  statement  It  reads,  It  takes 
ANOMLH  as  saro  in  this  case;  ANOMIA,  statements 
005^  and  0060,  on  the  other  hand,  is  operative, 
at  least  within  the  range  of  RATIO  >1.  In  a 
similar  manner,  the  exponential  factor  in  ANOMEE, 
statement  0Ct>6,  and  given  in  Eq.  (4)  , is  written 
in  statement  0062  as  CUTOFF.  The  electron-electron 
resistivity  is  turned  off  when  the  beam  current  den- 
sity CRDKNT  drops  below  the  minimum  value  required 
for  the  electron-electron  two  stream  instability, 
STABLJ  or  J . as  given  in  Eq.  (6) . This  same 
cutoff  also  i?mits  the  wave  resistivity,  ANOMWV, 

CSN  0067-006o,  although  in  this  run  it  is  set 
at  aero. 

CSN  006?-0099  Within  these  statements  are  set  the  time  increments 
for  which  the  computations  are  performed.  The 
calculation  time  step  length  HO,  is  set  at  2 nsec 
initially.  The  first  SO  time  steps  have  2 nsec 
Increments,  the  next  50  steps  have  5H0  or  10  nsec 
increments,  (CSN  00^0  and  00^4),  Additional 
steps,  if  more  than  100  were  called  for,  would 
have  increments  of  25HC  or  SO  nsec.  (CSN  00^3 
and  00Q6) . 

B.  Subroutine  FINI: 

TWO BOD  and  TRHEEB  are  the  two  and  three -body 
recombination  rate  coefficients  respectively. 

These  are  expressed  as  o3  and  ar.x  in  Eqs.  (1^) 
and  (l6)  in  Ref.  1. 

A NO ME  and  AN CM I are  the  electron  and  ion  portions 
of  the  resistivity  as  expressed  in  Eqs.  (12a)  and 
(12b),  Ref.  1,  respectively. 

The  equations  calculated  in  subroutine  FINI,  ENTRY 
F,  are  listed  in  statements  0062,  006?  and  0064, 
being  the  electron  energy,  electron  and  ion 
density,  and  ion  energy  equations  respectively. 

These  are  expressed  in  Eqs.  (^) , (10)  and  (11)  in 
Ref.  1. 

I NT,  EXTINT  and  ERROR  : 

I NT  chooses  an  Integration  accuracy  appropriate 
to  the  machine  on  which  it  is  computing,  or  1 
part  in  10s,  whichever  is  greater.  Then  it 
involves  subroutine  EXTINT,  as  necessary  to 
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CSN  00 '2,  00?? 


CSN  0057,  0058 


CSN  0062-0064 


C.  Subroutines, 
CSN  0005  (INT) 


T 


I 


CSN  0010  (INT) 
CSN  0018  (INT) 
CSN  0022  (INT) 


integrate  the  energy  and  density  equations  from  X to 
X + HO.  A description  of  the  arguments  of  INT 
appears  at  the  beginning  of  EXTINT.  COMMON /ERRC0M 
provides  subroutine  ERROR  with  variables  it  needs. 

Initialize  "S"  and  the  local  variables. 

Determine  an  appropriate  H. 

Perform  the  required  integration. 


D.  Computation  Print  Out: 

A sample  computation  printout  section  has  been  reproduced 
to  illustrate  the  existing  format  of  calculation  display.  Initial 
values,  as  previously  defined,  are  listed  across  the  top  of  each 
section,  and  beneath  these  are  the  15  columns  of  calculations. 

The  sample  calculation  has  initial  values  as  follows:  time  step 
2 nsec,  initial  electron  density  5 X 10i4cm_3,  initial  temperature 
2.5  eV,  maximum  beam  current  density  15  kA/cma,  neutral  fill 
density  1 x 1016cm*3  or  150  mT,  and  applied  magnetic  field  2 kG. 
Note  that  all  the  resistivities  are  operative  except  the  electron- 
electron.  Also  note  how  the  ion  acoustic  resistivity  is  "turned 
off"  at  about  75  nsec  when  RATIO  falls  below  0.3. 
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